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1. 


INTRODUCTION 


This  report  presents  a  numerical  method  for  computing 
the  effects  of  an  arbitrary  rotation  of  the  coordinate  system 
on  a  spherical  harmonic  expansion.  The  analysis  was  motivated 
by  a  study  of  the  approximations  involved  in  the  use  of  the 
Fourier  In-Flight  Transfer  Function  for  the  evaluation  of 
high-frequency  gravity  missile  impact  errors  (Ref.  6). 

The  method  has  a  wide  range  of  applications.  In 
Ref.  6  it  was  used  to  minimize  latitude  distortion  on  a  planar 
grid  representing  surface  gravity  in  the  neighborhood  of  the 
launch  point  by  redefining  the  coordinate  system  to  make  the 
(false)  equator  pass  through  the  launch  point.  Similar  appli¬ 
cations  exist  to  data  reduction  and  model  validation  where 
reference/model  gravity,  usually  described  by  a  spherical 
harmonic  expansion,  has  to  be  evaluated  on  a  local  grid  in  the 
region  where  data  exists.  Other  examples  of  possible  applica¬ 
tions  are: 

•  Generation  of  polar  charts  and  maps 
where  any  given  point  is  the  pole  of 
the  expansion 

•  Direct  simulation  and  gravity  model 
generation  for  false-pole  implementa¬ 
tions  of  inertial  navigation  and 
guidance  systems 

•  Comparison  and  correlation  studies  of 
different  geophysical  global  models 
(Ref.  3) 

•  Theoretical  analyses  of  isotropic  and 
stationary  properties  of  the  gravity 
field . 
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The  problem  of  determining  the  transformation  of  a 
spherical  harmonic  expansion  under  arbitrary  rotations  of  the 
coordinate  system  has  received  considerable  attention  in  the 
past.  Equivalent  formulas  have  been  independently  derived  by 
several  investigators.  The  earliest  reference  is  to  A.  Schmidt 
and  dates  back  to  1899  (Ref.  13).  Other  works  often  quoted  in 
the  literature  are  due  to  Herglotz  (Ref.  4),  Wigner  (Ref.  15) 
and  Jeffreys  (Ref.  9).  More  recently,  Kleusberg  (Ref.  12)  has 
derived  an  approximation  valid  for  small  rotations.  Applica¬ 
tions  to  Geodesy  and  the  study  of  satellite  orbits  have  been 
given  by  Kaula  (Ref.  10),  Cook  (Ref.  2),  Izsak  (Ref.  8), 

Baimino  and  Borderies  (Ref.  1)  and  Giacaglia  (Ref.  5). 

Because  of  the  complexity  of  the  general  transforma¬ 
tion  formula  and  because  of  numerical  instability  in  the  propa¬ 
gation  of  the  tranformation  coefficients,  the  formula  has  only 
been  used  for  low  degree  (<10)  expansions.  The  numerical 
method  presented  in  this  report  has  been  tested  and  yields 
accurate  results  with  expansions  to  degree  180.  The  method  is 
based  on  a  decomposition  of  an  arbitrary  rotation  of  the  co¬ 
ordinate  system,  R,  into  five  elementary  rotations  in  the  form 

R  =  A  B'1  Aq  B  A  (1-1) 

Y  p  u 

where  Afl ,  A^ ,  and  A^  represent  polar  axis  rotations  and  B  and 
B  ^  are  90-degree  rotations  about  a  fixed  equatorial  axis  whose 
effect  on  the  spherical  harmonics  is  described  by  a  set  of 
precomputed  weighting  coefficients. 

This  report  is  organized  as  follows:  Section  2  pre¬ 
sents  Lhe  derivation  of  Eq .  1-1  and  discusses  the  effect  of 

each  of  the  rotations  A  ,  A„ ,  A  ,  B  and  B  ^  on  the  spherical 

a  B  y 

harmonic  coefficients.  In  Section  3  an  algorithm  is  given  for 


the  evaluation  of  the  weighting  coefficients  associated  with 
rotations  B  and  B  ^ .  Section  A  presents  some  numerical  results 
obtained  with  expansions  to  degree  and  order  180.  Finally,  a 
summary  of  the  contents  of  this  report  is  given  in  Section  5. 
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2. 


ALGORITHM  FOR  ROTATION  OF  SPHERICAL  HARMONICS 


In  this  section,  a  discussion  of  two  different  con¬ 
ventional  definitions  of  spherical  harmonics  is  given  first. 
Simple  relations  are  given  for  transforming  coefficients  from 
one  convention  to  the  other.  This  discussion  is  followed  by  a 
presentation  of  the  general  formula  for  rotation  of  spherical 
harmonics.  Finally,  the  algorithm  for  practical  implementation 
of  this  transformation  is  introduced. 


Perhaps  the  most  widely  accepted  definition  of  normal¬ 
ized  surface  spherical  harmonics  is  that  of  Ref.  7: 


R£Q(e,A)  =  (2£+l  )1/2  P£0(cos6) 


for  £  >  0 ,  and 
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for  1  <  m  <  £,  where  P„  is  the  associated  Legendre  function 
—  —  £m 

of  degree  £  and  order  m  given  by 
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with  6  and  A  as  in  Fig.  2-1.  However,  the  equations  and  re¬ 
sults  in  this  report  are  considerably  simplified  when  complex 
spherical  harmonics  are  used.  The  normalized  complex  surface 
spherical  harmonics  are  defined  as  (Ref.  4) 


R  92994 
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Figure  2-1  Spherical  Coordinates 


yJ(6,\)  = 


<2i!+1>  rfeirf 


1/2  1W<:os  e> 


(2-3) 


where  i  =  V  -  1 . 

To  simplify  the  terminology,  in  the  sequel  the  func¬ 
tions  defined  by  Eq.  2-1  are  referred  as  "real  harmonics"  and 
those  given  by  Eq .  2-3  as  "complex  harmonics"  or  simply  as 
"harmonics"  when  the  meaning  is  clear  from  the  context. 

For  fixed  degree  £  and  nonzero  order  m,  there  are  two 

real  harmonics:  R„  and  Sn  .  This  is  also  the  case  for  the 

£m  £m  _m 

complex  harmonics.  The  functions  Y™  and  Y£  correspond  to 
order  ■  in  (  .  In  other  words,  the  definitions  in  Eqs .  2-2  and 
2-3  also  apply  to  negative  m  as  long  as  |m|<£.  It  can  be  shown 
that 


y£  (e  ,x )  =  (-l)“l  y“  (e  ) 


(2-4) 


where  the  superimposed  bar  denotes  complex  conjugation. 

The  expansion  of  a  function  h(6,A)  in  a  series  of 
complex  harmonics  is  of  the  form 


a»  £ 

h(6,A)  =  ^  Cj  Yj(e,A) 

£=0  m=-£ 


(2-5) 


where  the  coefficients  are  complex  numbers  given  by 


C£  =  J  h(0’A)  do 


(2-6) 


Since  all  functions  h(0,A)  considered  in  this  report  are  real, 
it  follows  from  Eqs.  2-4  and  2-6  that  the  coefficients  C^1  must 
satisfy 


C'm  =  (-l)m  c£ 


(2-7) 


Now,  suppose  that  the  expansion  of  h(0,A)  in  real  harmonics  is 


co  £ 


®  £ 


•»  =  S  E  “in,  IWe'X)  *  Z  Z  b*n.  S*m(6'X) 


£ =0  m=0 


£=1  m=l 


(2-8) 


Using  Eqs.  2-1,  2-3,  2-4,  and  2-5  it  can  be  shown  that 


C£  =  a£  0 


£  >  0 


(2-9) 


C£  =  <*£m  -  ihzJ^ 


1  <  m  <  £ 


(2-10) 


These  relations,  together  with  Eq .  2-4  permit  direct  transla¬ 
tion  of  expansions  in  real  harmonics  into  expansions  in  complex 
harmonics  and  vice  versa. 

Next,  consider  a  rotation,  R,  of  the  coordinate  system 
through  the  Euler  angles  a,  p,  y  defined  as  in  Fig.  2-2.  It 
is  shown  in  Ref.  4  that  a  harmonic  Y™(0,A)  is  mapped  into  a 
linear  combination  of  harmonics  of  the  same  degree  in  the  form 

£ 


Y®<e,A) 


£ 

j  =  -£ 


Q™,J  yJ(e  •  ,a-) 


(2-11) 


where  0'  and  A'  are  colatitude  and  longitude  in  the  rotated 
coordinate  system  and  where  the  weights,  depend  on  the 

Euler  angles.  Combining  the  formulas  for  the  weights  given  in 
Ref.  4  with  the  results  of  Szego  (Ref.  14)  on  extensions  of 
Jacobi  Polynomials,  P^vl*v2\  to  all  possible  values  of  the 
superindices,  it  can  be  shown  that 
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Figure  2-2  Euler  Angles.  The  angles  a,  p,  y  are  successive 
rotations  about  the  Z  axis,  the  nodal  axis,  and 
the  Z'  axis,  respectively. 


Using  Eq.  2-11,  the  expansion  in  the  original  coordi¬ 
nate  system,  Eq .  2-5,  can  be  transformed  into  an  expansion  in 
the  primed  system.  The  result  is 


w  JC 

EE 

£  =  0  m= - £ 


cj  Y®(e )  = 


-  A* 

EE 

£=0  j  =  -£  1_ 


E  - 

m=-£ 


m,  j 


Vj<6 


,A’  ) 


(2-13) 

The  expression  in  brackets  is  readily  recognized  as  the  coef¬ 
ficient  of  degree  £  and  order  j ,  ,  in  the  new  coordinate 

system.  Thus, 


£ 

t>i  =  T.  Q”lJ  c? 


(2-14) 


If  the  above  relation  is  viewed  as  a  computational 
equation,  the  transformation  weights,  Q^’^,  have  to  be  evalu¬ 
ated  each  time  a  new  choice  of  Euler  angles  is  used.  The  de¬ 
pendence  on  the  angles  a  and  y  presents  no  difficulty.  However, 
recursions  on  which  are  valid  for  a  range  of  values  of  p 

rapidly  become  numerically  unstable. 

The  alternative  suggested  in  this  report  is  to  decom¬ 
pose  the  general  rotation  of  Fig.  2-2  into  a  collection  of 
five  successive  rotations  each  of  which  has  transformation 
weights  which  are  either  dependent  upon  an  angle  but  easily 
computed  or  are  fixed  and  precomputed.  The  decomposition  and 
the  transformation  weights  associated  with  the  component  rota¬ 
tions  are  discussed  below. 

Next,  it  is  shown  that  the  general  rotation  of  Fig.  2-2 
can  be  decomposed  as 


R  =  A  B  A0  B  A 
Y  pa 


(2-15) 


where  A^  represents  a  rotation  by  an  angle  0  about  the  axis  of 
the  third  coordinate  (the  polar  axis),  and  B  is  a  rotation  of 
n/2  about  the  axis  of  the  second  coordinate.  (Figure  2-1  de¬ 
fines  the  ordering  of  the  axes.) 


The  appearance  of  the  first  and  last  terms  on  the 
right  side  of  Eq .  2-15  is  clear  from  Fig.  2-2.  They  corre¬ 
spond,  respectively,  to  the  rotations  about  the  transformed 
and  the  original  polar  axes  in  Fig.  2-2.  Consider  now  the 

factor  B  ^  A0  B  in  Eq .  2-15.  The  matrix  representations  of 
P 

A„  and  B  are 


and 


B  = 


0 

1 

0 


Therefore,  the  product  B_1A0B  is 

p 


f1  ° 

B  ^ApB  =10  cos(S 

\  0  -sinp 


(2-16) 


(2-17) 


(2-18) 


which  is  precisely  the  matrix  representation  of  the  intermediate 
rotation  in  Fig.  2-2.  This  shows  the  validity  of  the  decomposi¬ 
tion  formula  (Eq.  2-15). 


Thus,  the  algorithm  for  transformation  of  spherical 

harmonic  expansions  consists  of  successively  computing  the 

effect  of  rotations  A  .  B,  AD  ,  B-^  and  A  on  the  harmonic  coef- 

“  P  Y 

ficients.  The  transformation  weights  associated  with  each  of 
these  rotations  are  discussed  below.  Note  that  only  those 
coefficients  with  nonnegative  upper  index  need  to  be  computed 
at  each  stage  because  the  function  represented  by  the  expansion 
is  real  and,  therefore,  each  of  the  rotations  will  yield  coef¬ 
ficients  which  satisfy  Eq.  2-7. 
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The  transformation  weights  of  a  polar  rotation,  , 

$ 

are  a  simple  phase  change  applied  separately  to  each  of  the 
harmonic  coefficients.  If  and  D™  are  the  original  and 
transformed  coefficients  then 

Dm  _  eim$  cm  (2-19) 

This  can  be  demonstrated  directly  from  the  definition  of 
complex  harmonics,  Eq .  2-3,  by  redefining  the  longitude 
reference  as  A. '  =  \-<j>. 

The  transformation  weights  of  rotation  B  are  given  by 
Eq .  2-12  with  a  =  n/2 ,  0  =  n/2  and  y  =  -n/2.  The  expression 
for  these  weights,  which  are  denoted  by  T™’J  in  the  sequel, 
reduces  to 


J  j 

£ 


<-l)£'m  2’j 


p(m+j , j-m) 

4-j 


(0) 


(2-20) 


Thus,  all  t^,j  are  real.  These  weights  also  possess  useful 
symmetry  properties  which  are  given  below. 


Using  the  results  of  Ref.  14  it  can  be  shown  that  the 
Jacobi  polynomials  satisfy 


p(m+j . j -m) 
i-3 


(0) 


2J -m 


(m+j ,m- j ) 
r£  -m 


(0) 


(2-21) 


It  then  follows  from  the  last  two  equations  that 


1/2 


(2-22) 


Hence , 


T^’111  (2-23) 

which  implies  that  the  numerical  evaluation  of  T^’^  can  be 
limited  to  j>jn. 

Similarly,  it  can  be  shown  that 


T"m’j  =  (-l)£-j  tJ'J  (2-24) 

This  relation,  and  the  fact  that  T^,J  is  real,  reduce  the  com¬ 
putational  burden  of  the  effect  of  rotation  B  on  the  harmonic 
coefficients  by  a  factor  of  four.  As  before,  let  C™  and  Dj'  be 
the  original  and  transformed  coefficients.  Then,  from  Eqs.  2-14, 
2-23  and  2-24, 


£ 

oj  =  tJ-J  c°  *  ^'J  +<-ni'(j'fro>  cj) 

m=l 


(2-25) 


When  compared  with  Eq.  2-14,  the  above  relation  shows  a  factor 
of  two  reduction  in  the  number  of  multiplications  because  of 
the  disappearance  of  negative  indices  in  the  summation.  In 
addition,  since  the  quantity  in  brackets  is  either  real  or 
purely  imaginary,  another  factor  of  two  is  gained. 
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To  conclude  this  section,  consider  rotation  B 
Euler  angles  are  a  =  n/2 ,  p  =  -n/2,  and  y  =  -n/2.  It  then 
follows  from  Eq .  2-12  that  the  transformation  weights  of  B 
denoted  by  W^’J,  are  given  by 


w"1  ’  J 
W£ 


(21  \17 

■j  \ll+m) 

fe)v 


p(m+j,j-m) 


(2-26) 


Comparing  this  expression  with  Eq .  2-22,  it  can  be  seen  that 
the  weights  of  B  are  identical  to  the  weights  of  B  except 
for  transposition  of  the  upper  indices;  i.e., 


»  j  —  'I' j  »  ® 

£  ‘  £ 


(2-27) 


Therefore,  a  single  table  of  weights  is  sufficient  to  describe 
the  effect  of  both  B  and  B  1  on  the  harmonic  coefficients.  In 
terms  of  T™’-*,  the  effect  of  B  *  is 


=  T^’0  C?  +  ^  Ti’m  [cT  ♦  (-l)£"(j+m)  C®]  (2-28) 


% 


3. 


EVALUATION  OF  THE  WEIGHTS 


The  transformation  weights,  T^'^,  are  evaluated  through 
recursion  on  the  three  indices.  All  recursion  relations  used 
in  the  algorithm  given  below  were  derived  from  formulas  in 
Ref.  14  relating  contiguous  Jacobi  polynomials.  The  slowest 
varying  index  in  the  algorithm  is  £.  Thus,  all  weights  T™’J 
for  fixed  degree  £  can  be  stored  as  a  single  record.  Later, 
when  the  table  of  weights  is  used  to  rotate  a  spherical  har¬ 
monic  expansion,  these  weights  can  also  be  read  as  a  single 
record  and  all  the  coefficients  of  degree  £  can  be  processed 
at  once. 


The  weights  of  degree  0  and  1  serve  as  starting  values 
for  the  computation.  They  are  given  in  Table  3-1  which  also 
lists  all  the  weights  through  degree  five.  Only  the  weights 
T£,J  for  0  <_  m  £  j  <_  £  need  to  be  computed  because  of  the  sym¬ 
metry  of  the  weights  (Eqs.  2-23  and  2-24)  and  because  the  algo¬ 
rithm  for  rotation  of  spherical  harmonic  coefficients  only  has 
to  evaluate  the  transformed  coefficients,  ,  for  j  0 . 

The  following  notation  is  used:  Int(x)  denotes  the 
largest  integer  less  than  or  equal  to  x  and  is  defined  as 

=  [ (£+k+l)(£-k) ]1/2  (k=0 . £ )  (3-1) 

Next,  the  algorithm  is  presented.  A  discussion  of  its  various 
steps  is  given  at  the  end  of  this  section. 
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TABLE  3-1 

WEIGHTS  T^1  ’ ^  FOR  £<5 


E 

m 

— 

j=o 

j=l 

■a 

B 

j=4 

■B 

1 

— — — — — 

— 

D 

0 

-72/2 

B 

73/2 

1/2 

2 

-1/2 

0 

76/4 

2 

0 

-1/2 

-1/2 

2 

2 

76/4 

1/2 

1/4 

3 

0 

75/4 

0 

-75/4 

3 

-75/4 

-1/8 

7To/8 

775/8 

3 

2 

0 

-710/8 

-1/2 

-76/8 

3 

3 

75/4 

715/8 

76/8 

1/8 

4 

0 

3/8 

0 

- 770/8 

0 

775/16 

4 

1 

0 

3/8 

72/8 

-77/8 

-774/8 

4 

? 

•m. 

-710/8 

-72/8 

1/4 

774/8 

77/8 

4 

3 

0 

-77/8 

-774/8 

-3/8 

-72/8 

4 

4 

770/16 

774/8 

77/8 

72/8 

1/16 

5 

0 

0 

-730/16 

0 

735/16 

0 

-377/16 

5 

1 

730/16 

1/16 

-77/8 

-742/32 

727/16 

7210/32 

5 

2 

0 

77/8 

1/4 

-73/16 

-73/4 

-730/16 

5 

3 

-735/16 

-742/32 

76/16 

13/32 

972/32 

375/32 

.. 

4 

0 

-771/16 

-73/4 

-972/32 

-1/4 

-770/32 

5 

— i 

lj 

377/16 

7210/32 

730/16 

375/32 

770/32 

1/32 
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Algorithm  for  Evaluation  of  the  Weights  T*‘1,J  (0  £  m  £  j  <  £  <  L) 


(Initialization] 


Let  £  =  2,  Tq,0=1,  T°’°=0,  T®,1=-f2/2  and  T^,x=  1/2 


1,1. 


In  addition  set  a«-T^’^,  and  d*--T^’^ 


[Determine  and  T^’^] 


If  £  is  even  then  set  a  «-  -  a  an<^  Put  =  a»  1 


If  9  is  odd  then  set  b  <- 
T®'1=b 


£(£-2 

U+1)U 


y 


b  and  put  t£’°=0. 


3.  [Determine  T£ for  j=2. 


£-1] 


Use  the  following  formula  recursively 


n0,j  _ 


A _  t°  » j ” 2 

J-l  l9 


(3-2) 


4.  [Determine  and  ^] 


M 


c  and  d 


,  _  _0 ,£  „1 ,£-l  , 

and  put  T.  =  c ,  T.  =  d 


Xs  O 


5.  (Determine  T^’£  for  m=l ,  ....  £ ] 


Use  the  follov-'n^  forme.1  a  refer  .  t .  « 


Tm  ,£  „ 

£ 


6.  (Determine  ^  for  m=2,  £-1] 


Use  the  following  formula  recursively 


Tm , £ - 1  _ 
SL 


m  SL  -m+ll  -jn-l.i-l 

m-1  £+m  j  l£ 

»  -J 


SL  -m+1 
£+m 


—m-1  ,£ 
SL 


(3-3) 


7.  (Determine  T^,J  for  2£m£j££-2] 
7.1  Set  k  =  max  { 0 ,£ -2- Int ( 2 (£-1 )1/2 ] } 


7.2  (Case  j-m  k] 


For  each  n  =  £-3,  £-4,  ...»  k  let  j  =  £-2,  £-3,  ....  n+1 
and  set  m  *■  j-n.  Compute 


Tm*  j 
£ 


m  m-1 

m  m£  _m+l,j  +  l  m  M£  -m-l.j+l 

3+i  „J  *  '  j+i  * 

m£  m£ 
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7.3  [Case  j-m<k] 


For  each  n=k-l ,  k-2  0  set  q  *■  Int  {  [  n+(  2£  ^-n^ )  1/2 } 

and  perform  the  following  two  steps  in  the  order  given 

7.3.1  For  j=q+l,  q+2 ,  ...»  £-2  set  m<-j-n  and  compute 

m-2  j 

j  -m+1  ^£ -an- 2  ,  j  _  2j  _  £  ^jm-l.j+l 

j+m-1  m-1  x£  j+m-1  m-1  £ 

M£  ^£ 

(3-6) 

7.3.2  For  j=q,  q-l»  •••,  n+1  set  m<-j-n  and  compute 

using  Eq.  3-5 

8.  [Update  degree] 

All  weights  of  degree  i  have  been  computed.  If  £<L  then 
set  £<-£  +  1  and  go  to  Step  2. 

Figure  3-1  presents  a  graphical  interpretation  of  the 
operation  of  the  algorithm.  The  numbers  shown  in  parentheses 
next  to  the  arrows  are  the  steps  of  the  algorithm  in  which  the 
weights  along  the  lines  parallel  to  the  arrows  are  computed. 

The  direction  of  the  arrows  indicates  the  progression  in  which 
the  weights  are  evaluated. 

Four  quantities  (denoted  by  a,  b,  c,  and  d  in  the 
algorithm)  are  needed  to  propagate  the  weights  from  degree  £-1 
to  degree  £.  They  are  initialized  in  Step  1  and  updated  in 
Steps  2  and  4  where  they  are  also  used  to  evaluate  four  weights 
of  degree  £:  T°’°,  T®’1,  T°,£  and  T*’*"1.  These  weights, 

indicated  by  circles  in  Fig.  3-1,  serve  as  "seeds"  for  the 
computation  of  the  remaining  weights  of  degree  £. 


•pro  y  j 
£ 
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Figure  3-1  Computation  of  the  Weights  of  Degree  £ 

The  main  recursion  formula  is  Eq.  3-5.  It  permits 
the  propagation  of  the  weights  along  lines  parallel  to  the 
diagonal  (Fig.  3-1).  Its  application,  however,  requires  prior 
evaluation  of  the  weights  along  the  boundary  lines  m=0,  j=£-l, 
and  j=£.  This  is  achieved  in  Steps  3,  5  and  6  of  the  algorithm. 

The  main  recursion  formula,  Eq .  3-5,  is  numerically 

2  2  2 

unstable  outside  the  circle  m  +j  £  (See  Fig.  3-1).  This 
i s  the  reason  for  dividing  Step  7  into  two  separate  cases.  If 
the  weights  to  be  computed  along  the  subdiagonal  j -m=n (=constant ) 
ai]  lie  inside  this  circle  then  Eq.  3-5  is  used.  If  they  do 
not,  then  another  recursion  formula,  Eq .  3-6,  is  necessary. 

The  critical  value  of  n  is  computed  in  Step  7.1  and  is  denoted 
bv  k  ir  the  algorithm.  Thus,  in  the  region  j-m^k  the  main  re¬ 
cursion  formula  is  used  (Step  7.2)  while  for  j-m<k,  Eq .  ^-t  is 


used  outside  the  circle  (Step  7.3.1)  and  the  main  recursion 
formula  is  used  inside  the  circle  (Step  7.3.2). 


Double  precision  (53  bit  mantissa)  was  used  in  an  IBM 
4341  computer  to  generate  all  weights  to  degree  180.  CPU  time 
was  approximately  3  minutes.  The  accuracy  of  the  computed 
weights  was  later  verified  using  recursion  formulas  independent 
of  those  used  in  the  algorithm  described  above.  These  formulas 
are 


-pi  » j 
£ 


T0  ,  j 
£ 


if  £  - j  is  even 


T0,  j  +  1 
£ 


if  £  - j  is  odd 


(3-7) 


and 


.-in.j  o  (£-m-l)  (£-j  - 1  )  +  (£ -m-j  - 1 )  (j+m+3  )  Tm+l,j+l 

*£  Z  j+m+3  j  m  £ 

m£  h£ 

j+1  m+1 

j+m+1  m£  M£  Tm+2,j+2  (3-8) 

j+m+3  j  m  £ 
h£h£ 


Equation  3-7  permits  the  computation  of  all  weights 
on  the  line  m=l  in  terms  of  those  for  m=0.  Note  that  T£,J  is 
the  last  weight  evaluated  along  each  subdiagonal  using  the 
main  recursion  formula  in  Steps  7.2  and  7.3.2  (See  Fig.  3-1). 
Therefore,  a  comparison  of  the  weights  T^’^  generated  by  the 
algorithm  and  those  computed  on  the  basis  of  Eq .  3-7  provides 
a  measure  of  the  numerical  behavior  of  the  main  recursion 
formula . 


T 


On  the  other  hand,  Eq .  3-8  gives  T^1  ’  ^  in  terms  of  two 
neighbors  on  the  same  subdiagonal.  This  formula  is  not  numer¬ 
ically  stable.  T  t  '  'as  used  on1  v  -■  .  •  «  -  V  on  numerical  closure 

by  repeatedly  substituting  weights  generated  by  the  algorithm 
for  >  J +  ^  ancj  jJ-n+2,j+2  Qn  t^e  right  hand  side  of  Eq .  3-8  and 

comparing  the  result  with  the  value  of  T^>'1  produced  by  the 
algori thm . 

The  tests  showed  that  the  accuracy  of  the  computed 
weights  degraded  with  decreasing  absolute  magnitude  of  the 
weights  themselves.  No  comparison  showed  less  than  six  (deci¬ 
mal)  significant  figures  of  agreement.  Since  the  order  of  mag¬ 
nitude  of  the  weights  ranges  over  five  decades  and  since  the 
smaller  the  magnitude  of  the  weight,  the  smaller  its  effect  in 
the  evaluation  of  harmonic  coefficients  (see  Eq .  2-25),  the 
results  obtained  with  the  algorithm  were  considered  more  than 
adequate.  Further  proof  of  the  validity  of  the  computation  of 
the  weights  is  presented  in  the  next  section. 
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4. 


NUMERICAL  RESULTS 


The  algorithm  for  rotation  of  spherical  harmonics  was 
tested  using  three  synthetic  expansions,  all  to  degree  180. 

In  each  expansion,  the  coefficients  were  zero-mean  normal  vari¬ 
ates  chosen  according  to  Kaula's  rule  for  the  degree  variances 
of  the  geopotential  (Ref.  10).  The  population  variance  for 
each  of  the  coefficients  of  degree  £  in  the  expansions  in  real 
spherical  harmonics  was 

o£  =  10"10/£4  £  >  0  (4-1) 

The  series  were  first  rotated  through  random  Euler 
angles  and  then  transformed  back  to  the  original  coordinate 
system.  The  resulting  and  the  original  coefficients  were  com¬ 
pared  and  differences  were  determined.  Double  precision  was 
used  throughout  the  computation. 

Figures  4-1  and  4-2  summarize  the  results.  In  Fig.  4-1 
a  plot  is  given  of  the  RMS  error  as  a  fraction  of  the  expected 
value  of  the  magnitude  of  the  coefficients.  The  plot  was  ob¬ 
tained  by  combining  the  results  of  the  three  test  cases.  Thus, 
in  evaluating  the  RMS  error  in  the  coefficients  of  degree  £, 

6£  +  3  differences  were  considered.  The  resulting  value  was  nor¬ 
malized  by  the  expected  magnitude  of  the  coefficients  given  by 


k 


£ 


2o£//2n 


(4-2) 


with  as  in  Eq .  4-1.  The  purpose  of  this  normalization  was 
to  compensate  for  the  decreasing  magnitude  of  the  coefficients 
in  order  to  provide  a  measure  of  the  number  of  exact  signifi¬ 
cant  figures  which  can  be  expected  from  the  application  of  the 
algorithm.  Thus,  Fig.  4-1,  indicates  that  up  to  degree  180 
the  standard  deviation  of  the  errors  is  always  more  than  thir¬ 
teen  orders  of  magnitude  smaller  than  the  average  size  of  the 
coefficients . 

Figure  4-2,  on  the  other  hand,  is  a  plot  of  worst- 
case  results.  It  shows  the  largest  relative  error  found  in 
all  three  test  cases  as  a  function  of  degree.  In  other  words, 
for  given  degree  £  ,  the  largest  ratio  between  the  absolute 
values  of  the  error  and  the  true  (original)  coefficient  in  all 
test  cases  was  plotted  in  Fig.  4-2.  It  was  observed  that  the 
•largest  relative  errors  tend  to  occur  in  those  coefficients 
whose  magnitude  is  small  compared  with  other  coefficients  of 
the  same  degree.  This  is  because,  in  general,  all  of  the  co¬ 
efficients  of  degree  £  enter  in  the  computation  of  any  single 
coefficient  of  the  same  degree  in  the  rotated  expansion. 
Therefore,  the  absolute  errors  have  the  same  order  of  magni¬ 
tude  for  all  computed  coefficients  but  they  represent  a  larger 
percent  in  the  smaller  coefficients.  As  Fig.  4-2  indicates, 
in  no  case  did  the  computed  coefficients  have  less  than  nine 
significant  figures  of  agreement  with  the  original  coefficients. 

3 

The  number  of  operations  in  the  algorithm  is  0(N  ) 
where  N  is  the  highest  degree  in  the  expansion.  In  an  IBM  4341 
computer,  CPU  time  for  a  single  (one-way)  rotation  was  slightly 
more  than  two  minutes  for  an  expansion  to  degree  180  and  less 
than  three  seconds  for  an  expansion  to  degree  36. 
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5. 


SUMMARY  CONCLUSIONS  AND  RECOMMENDATIONS 


An  algorithm  has  been  given  for  accurately  computing 
the  effect  of  an  arbitrary  rotation  of  the  coordinate  system 
on  a  spherical  harmonic  expansion.  The  method  consists  of 
decomposing  the  rotation,  R,  into  5  successive  component  rota¬ 
tions  in  the  form 

R  =  ay  B-1  Ap  B  Aq  (5-1) 

where  or,  p  and  y  are  the  Euler  angles  of  R. 

The  effect  of  A  ,  Aft  and  A  on  the  harmonic  coeffi- 
cients  is  easily  computed  since  they  represent  rotations  about 
the  polar  axis.  On  the  other  hand,  matrices  B  and  B_1  corre¬ 
spond  to  a  90  degree  rotation  about  an  equatorial  axis  and  its 
inverse.  Their  effect  on  the  harmonic  coefficients  was  shown 
to  be  described  by  a  set  of  weights  which  is  fixed  and  inde¬ 
pendent  of  rotation  R.  A  stable  algorithm  was  given  for  the 
evaluation  of  the  weights. 

Numerical  tests  with  synthetic  expansions  to  degree 
180  were  performed.  These  tests  verified  the  methodology  and 
demonstrated  the  high  accuracy  of  the  algorithm. 

It  is  recommended  that  the  algorithm  be  implemented 
in  the  Weapons  Support  System  (WSS)  as  a  general  tool  for  the 
analysis  and  use  of  spherical  harmonic  geopotential  models. 
Specific  and  immediate  applications  include  data  reduction, 
model  validation,  charting  and  false-pole  gravity  model 
generation . 
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